In [24]:
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
from scipy.integrate import odeint
from IPython.display import Image
from IPython.core.display import HTML
%matplotlib inline
A possible method to obtain this relationship is to integrate the reaction rate, $r_{\alpha}$, for the impurity production over time. \begin{equation} \alpha = \int_0^t \frac{d \alpha}{dt}dt \equiv \int_0^t r_{\alpha}(\{p_i\}_{i=1}^N; T, H)dt \end{equation} Here, $\{p\}_{i=1}^N$ denotes a set of N other possible parameters that the reaction depends upon (which will be denoted as $\{p_i\}$ for short). Given a model for the reaction rate, $r_{\alpha}$, the above equation can be integrated to give the impurity as a function of time, which can be used to determine an approximate shelf life. The difficulty is thus estimating an accurate kinetic model. Fortunately, observations of many reactions have lead to several assumptions help simplify the estimation of a kinetic model.
One assumption is that the temperature dependence of the reaction kinetics can be separated from the concentration dependence, and takes the form of an exponential function, which was first proposed by Arrhenius. The ASAP literature further suggests thatthe humidity dependence can also be separated as an exponential function. Thus, the temperature and humidity are often separted out into a rate constant, $k$. \begin{equation} r_{\alpha}=k(T,H;A,E)f(\{p_i\}) = A \exp^{-\frac{E}{RT}+BH}f(\{p_i\}) \end{equation}
The parameters, $p_i$, are often other reactants, which can be similarly expressed using reaction kinetic models. Therefore, the solution is obtained by the integration of the coupled reactions for all the species in the reaction. \begin{equation} \frac{d p_i}{dt} = k_i \hat{p_i}(\{p_j\}) \end{equation} Since the tolerable level of impurity is usually low, it is reasonable to assume that the concentration of the components involved in the reaction do not deviate substantially from the initial values. Thus, the above reactions can be approximated by a truncated Taylor series. \begin{equation} \hat{p_i} \approx \sum_j \frac{\partial p_i}{\partial p_j}\bigg|_{j \ne i}(p_j-p_{j,o}) = \sum_j \nu_{ij}(p_j-p_{j,o}) \end{equation} Here, $\nu_{ij}$ is the stoichiometric coefficient ratio between reactants $i$ and $j$. Substiting this into the the coupled equations. \begin{eqnarray} \frac{d \alpha}{dt} = k_\alpha \nu_{\alpha j} \sum_j (p_j-p_{j,o}) \\ \frac{d p_i}{dt} = k_i \nu_{ij} \sum_j (p_j-p_{j,o}) \end{eqnarray}
Integration of the above equations appxroximates the changes of impurity concentration for small step changes, and is exact if all the kinetics are first order. Moreover, the equations can be easily solved numerically or can be analytically solved using an orthogonal tranformation. Thus, this is a reasonable starting point for reactions with sufficient measurement data, such that the reaction mechanism (i.e. the $\nu_ij$) can be determined, as well as the temperature and humidity dependence parameters ($\{E_i\}$, and $\{B_i\}$) be estimated. However, this may not be suitable for the impurity analysis, which has very limited data.
The idea of isoconversion is to esitmate the temperature and humidity parameters without developing a mechanistic model. First, the imprutiy reaction rate is written by separating out the temperature and humidity from the other reactant dependencies. \begin{equation} \frac{d \alpha}{dt} = A \exp^{-\frac{E}{RT}+BH}f(\{p_i\}) \end{equation} Next assume that each reactant species can be parametized by the impurity concentration. \begin{equation} p_j = \hat{p_j}(\alpha) \end{equation} An example of a parameterization often used in kinetics courses is that the concentration of the reactant decreases proportionally (stoichiometrically) to the product (i.e. $p_j = p_{j,o} - \nu_{j\alpha} \alpha$). Substituting the general parametrized relationships into the reaction equation gives. \begin{equation} \frac{d \alpha}{dt} = A \exp^{-\frac{E}{RT}+BH}f(\{\hat{p_j(\alpha)}\}) \equiv A \exp^{-\frac{E}{RT}+BH} f_\alpha(\alpha) \end{equation} If a series of experiments are conducted at different temperature and humidity trajectories, but the instantaneous reaction rate, $d\alpha/dt$, is recorded at identical impurity concentrations, $\alpha_{sp}$, then the data can be regressed to obtain the modified Arrhenius parameters. \begin{equation} \ln \left( \frac{d \alpha}{dt} \right )_{\alpha_{sp}} = \ln A - \frac{E}{RT} + BH + \ln f_\alpha(\alpha_{sp}) = \ln(A f_\alpha(\alpha_{sp})) - \frac{E}{RT} + BH \end{equation} As discussed in literature, estimating the parameters using data for the instantaneous rate at level $\alpha$ requires differential type measurements, and is such suited for measurements using DSC and similar. For experiments that measure the net change of impurity level (such as TGA, or purity tests after set durations) the integral form of the isoconversion analysis is more suitable.
Integrating the general relationship between rate and parameters gives. \begin{equation} \int_0^\alpha \frac{d \alpha}{f_\alpha(\alpha)} \equiv g(\alpha) = \int_0^t A \exp^{-\frac{E}{RT}+BH} dt = A \exp^{-\frac{E}{RT}+BH}t \end{equation} Finally, rearranging renders the desired result. \begin{equation} \ln t = \frac{E}{RT} + B H + \ln \left[\frac{g(\alpha)}{A}\right] \end{equation} The above gives the expression for time, $t$, to reach an impurity concentration, $\alpha$ for a specified temperature, $T$, and humidity, $H$. The idea is to make a series of measurements at different conditions ($T$, $H$) holding the samples until they achieve equal degredation, such that $g(\alpha)$ remains constant, and thus leave only $E$, $B$, and $g(\alpha)/A$ as unknowns. Once the three parameters are regressed from experimental data, the shelf life can be estimated.
The assumptions of the isoconversion method are:
It is unlikely that intermediate and competing reactions will not be temperature dependent. However, the second assumption can be modified to allow some temperature and humidity dependence. Assume that the dependence of the degredation rate on the other reactants is completely multiplicative. \begin{equation} f(\{p_i\}) = \Pi_i k_i(T,H) p_i \end{equation} where $k_i$ is of the form of the modified Arrehnius. If the $p_i$ can be parametrized by $\alpha$, then the degredation reaction rate will still be of the form. \begin{equation} \frac{d \alpha}{dt} = k^*(T,H) \Pi_i \hat{p_i} = A^* \exp^{-\frac{E^*}{RT}+B^*H} f_\alpha^*(\alpha) \end{equation} Here the asterics denote an effective value from all the reactions, but is inconsequential in the analysis. Lastly, if the amount of degredation is low, than the functional dependence of the degredation reaction on other species can be expanded in a Taylor series. \begin{equation} f(\{p_i\}) \approx \sum_j k_j(T,H) \frac{\partial f(\{p_i\})}{\partial p_j} \bigg |_{j \ne i}(p_j - p_{j,o}) \end{equation} Again, assuming the rate constants, $k_j(T,H)$, are a modified Arrhenius form and the $p_j$ can be parametrized by $\alpha$. \begin{equation} \frac{d \alpha}{dt} = A \exp^{-\frac{E}{RT}+BH} \sum_i k_i(T,H) f_\alpha(\alpha) \end{equation} The above can be simplified if all the $k_i(T,H)$ are identical (not likely) or one of the rate constants is much greater than the other (more likely). In this scenario, the equation reduces to the same form as in the previous assumptions. \begin{equation} \frac{d \alpha}{dt} = A^* \exp^{-\frac{E^*}{RT}+B^*H} f_\alpha^*(\alpha) \end{equation}
Lastly, if the the time dependence of the individual $p_i$ are coupled as a series of ordinary differential equations(as was demonstrated in the mechanicstic section). \begin{equation} \frac{d p_i}{dt} = k_i \hat{p_i}(\{p_j\}) \end{equation} If one the equations is rate limiting and the dependence of the reactants in this equation can be parametrized by the degredation extent, then the result reduces to the same conclusion as seen in the previous two assumptions.
This added analysis shows that the previous assumption 2 can be augmented to include temperature dependence for degredation reactions that have a multiplicative dependence on other reactants, for small degredation extents with a dominant pathway, and for a series reaction where one intermediate is rate limiting. One of these assumptions would likely be satistifed degredation reactions dominated by a particular pathway over all temperature and humidity enviroments measured. Contrarily, the assumptions would be invalid for competing reactions that have strong temperature dependence over the range of experiments.
The following examples are executable calculations that can be used to evaluate various models.
In this example the impurity product, $P$, can either be produced from crystalline drug, $D$, or an amorphous intermediate, $I$, that is more reactive (this example is extracted from Waterman et al. Pharm Res Vol 24 2007 p 783).
Expressing the system in matrix form. \begin{equation} \frac{d}{dt} \begin{bmatrix} P \\ I \\ D \end{bmatrix} = \begin{bmatrix} 0 & k_2 & k_3 \\ 0 & (-k_{1r} - k_2) & k_{1f} \\ 0 & k_{1r} & (-k_{1f}-k_3) \end{bmatrix} \begin{bmatrix} P \\ I \\ D \end{bmatrix} \end{equation}
Below is a numerical integration of the the above coupled ODEs. The kinetic parameters and initial conditions can be changed, and the plot will update accordingly. To edit the parameters:
In [62]:
"""Don't edit anything here unless you know what you're doing."""
# time domain
total_time = 1000
time_increment = 1000
t = np.linspace(0, total_time, time_increment)
# definition of the derivative of the vector variable.
def deriv(y, t, k1f, k1r, k2, k3):
P, I, D = y
dPdt = k2*I + k3*D
dIdt = (-k1r-k2)*I + k1f*D
dDdt = k1r*I + (-k1f-k3)*D
dydt = [dPdt, dIdt, dDdt]
return dydt
In [75]:
1e10*np.exp(-20/(0.00198588*298))
Out[75]:
In [71]:
"""User defined Arrhenius constants and initial conditions.
See the reaction diagram directly above for variable definitions."""
# k1
A1f = 1e10
E1f = 20
k1f = 1e-4
k1r = 1e-5
k2 = 1e-4
k3 = 1e-3
Po = 0
Io = .1
Do = 0.9
yo = [Po, Io, Do]
# call solver
yt = odeint(deriv, yo, t, args=(k1f, k1r, k2, k3))
plt.plot(yt)
ax = plt.gca()
ax.set_xlabel("time (s)")
ax.set_ylabel("concentration (%)")
ax.legend(["P", "I", "D"])
Out[71]:
In [70]:
# calculate impurity totals from different contributions
P = yt[:,0]
I = yt[:,1]
D = yt[:,2]
r = D*k3
PD = np.array([sum(r[0:i]) for i in xrange(len(r))])
PI = P - PD
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t, PD)
ax.plot(t, PI)
ax.plot(t, P)
ax.legend(["impurity from D", "impurity from I", "total impurity"], loc="best")
Out[70]: